#ifndef CURRENT_H
#define CURRENT_H
#include "cdecl.h"
/**
 * The module implements current abstraction. 
 * In our calculations we have I(t,l) - in Ampers for every trajectory 
 * point at any time moment. A time spectrum I(w,l) of this function is 
 * present in numeric calculations.
 * 
 * We define signal spectrum this way:
 * S(w) = \frac{1}{2\pi}\int\limits_{-\infty}^{+\infty} S(t) e^{+iwt} dt
 */

/**
 * Our currents are like this:
 * 
 * I(l,t) = \Phi(t-\tau(l)) \ksi(l)
 * 
 * Where \tau is something like l/v and ksi - something like exp(-eta*t).
 * If and only if ksi is a constant, there are no charges left after 
 * impulse on the trajectory. 
 * 
 * \tau function is assumed to be monotonically increasing (insane otherwise)
 * \tau(0) = 0
 * 
 * \ksi function is monotonocally decreasing 
 * \ksi(0) = 1 
 */
enum current_shape { // The \Phi function
	CURRENT_SHAPE_EXPEXP, // Two exponents: \Phi(\tau) = I0 (exp(-\alpha*\tau) - exp(-\beta*\tau))
	CURRENT_SHAPE_POWEXP  // Power of t and an exponent: \Phi(\tau) = I0 \tau^{n} exp(-\alpha*\tau)
};

struct current_shape_expexp {
	enum current_shape id;
	double Q; // Integral over impulse - the carried charge
	double alpha; // 
	double beta; // beta > alpha, otherwise current becomes negative
};

struct current_shape_powexp {
	enum current_shape id;
	double Q;
	double alpha;
	uint8_t n; // n > 0, n mustn't be too great
};

union current_shape_u {
	enum current_shape id;
	struct current_shape_expexp expexp;
	struct current_shape_powexp powexp;
};




enum current_delay { // The \tau part
	CURRENT_DELAY_LINEAR, // The simple Dalamber-way t = l/v
	CURRENT_DELAY_EXPONENTIAL // v = v0 exp(-gamma * t) -> t = -1/gamma*ln(1 - l * gamma/v0) // gamma > 0, gamma = 0 is a singularity
};

struct current_delay_linear {
	enum current_delay id;
	double v0;
};

struct current_delay_exponential {
	enum current_delay id;
	double gamma;
	double v0;
};

union current_delay_u {
	enum current_delay id;
	struct current_delay_linear linear;
	struct current_delay_exponential exponential;
};




enum current_decrease { // The ksi part
	CURRENT_DECREASE_NONE, // ksi = 1
	CURRENT_DECREASE_EXPONENTIAL // ksi = exp(-eta * l)
};

struct current_decrease_none {
	enum current_decrease id;
};

struct current_decrease_exponential {
	enum current_decrease id;
	double eta;
};

union current_decrease_u {
	enum current_decrease id;
	struct current_decrease_none none;
	struct current_decrease_exponential exponential;
};


struct current {
	union current_shape_u		shape; // \phi
	union current_delay_u		delay; // \tau
	union current_decrease_u	decrease; // \ksi
};

/**
 * Calculates current spectrum I_w(l) for given frequency and distance
 * along discharge path.
 * @param in_current - current structure
 * @param w - 2*pi*frequency
 * @param l - length from discharge root
 * @return spectrum value
 */
complex current_get_spectrum(const struct current * in_current, double w, double l);

/**
 * For all lengths in [0,l] finds T such that most of the impulse energy
 * fits in [0,T]. This is needed to find such times that there is appropriately
 * no signal energy outside the diapason for the whole trejaectory.
 * @param cur - current object to study
 * @param l - trajectory length
 * @param eps - cut-off criteria: if energy_remains / full energy < eps, that's
 * enough to believe that the whole signal fits in [0,T]
 */
double current_get_duration(const struct current * cur, double l, double eps);

/**
 * Checks if the parameters hold correct values.
 * @param cur - current to check to be corretc
 * @return pointer to string with user-friendly description of problem,
 * that's does not require free-ing it, or NULL if all OK. 
 */
const char * current_validate(const struct current * cur);

/**
 * Finds such frequency that overall current spectrum for that frequency
 * is less than eps related to maximum value.
 * 
 * Generally speaking, we shold have integrated over the requencies here... But it's tricky. May be, numerically.
 * 
 * @param cur - current to stduy
 * @param eps - precision
 * @return the w frequency border (2*pi*f, [Rad/s])
 */
double current_get_max_w(const struct current * cur, double eps);


/**
 * By given current current object creates another one such that 
 * I_refl(l,t) = (-1) * I(l,t)
 * @param src - initial current object to reflect
 * @param refl - reflected current
 */
void current_mirror(const struct current * src, struct current * refl);

#endif
